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ABSTRACT 


EFFICIENT GRID GENERATION 
by 

Bill Gutierrez 

Master of Science in Physics 

Because the governing equations in fluid dynamics contain partial differentials 
and are too difficult in most cases to solve analytically, these differentials axe 
generally replaced by finite difference terms. These terms contain terms in the 
solution at nearby states. This procedure discretizes the field into a finite number 
of states. These states, when plotted, form a grid, or mesh, of points. It is at 
these states, or field points, that the solution is found. 

The optimum choice of states, the x,y,z coordinate values, minimizes error and 
computational time. But the process of finding these states is made more difficult 
by complex boundaries, and by the need to control step size differences between 
the states, that is, the need to control the spacing of field points. 

One solution technique uses a different set of state variables, which define a 
different coordinate system, to generate the grid more easily. A new method, de- 
veloped by Dr. Joseph Steger, combines elliptic and hyperbolic partied differential 
equations into a mapping function between the the physical and computational 
coordinate systems. This system of equations offers more control than either equa- 

viii 


tion provides alone. The author has modified Steger’s algorithm in order to allow 
bodies with stronger concavities to be used, offering the possibility of generating 
a single grid about multiple bodies. Work has also been done on identifying ar- 
eas where grid breakdown occurs. This work was supported by NASA under the 
NASA Minority Graduate Researcher’s program. 
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1 Introduction 

The governing equations in fluid dynamics, as in many other fields, contain terms 
in partial differentials. In the numerical solution of these equations the differen- 
tials are generally replaced by finite difference terms derived from taylor series 
expansions. The solution is then found from given auxiliary data. For instance, 
if u is a function of n variables 

u = u(x i,aJ 2 ,...,*„) (1) 

and the equation relating u to these variables is a partial differential equation, 
each m’th partial derivative of u can be expressed as a sum of terms in u, and an 
error term 

d m U 1 

^jT — ^ , . . . X{ “I - bj /A X i , . . . Xu) O ( (2) 

ax i j = 1 

where a,j and bj are constants and 0(Ax) p is the error term. It is important to 
note that, similar to the exact differential on the left, all other x variables are 
held constant while x; varies to produce the different u values. 0(Ax) p , the order 
of error, is a measure of difference between the exact differential and the finite 
difference term. The smaller Ax is, the smaller the error. Because Ax is generally 
less than one, higher values of p will also decrease this error. 

A more specific example of the foregoing is a finite difference equation for the 
first partial derivative of u at x t - = x\. Dropping notation of x variables held 
constant, equation (2) becomes in this case 


g W ) , .„(*') = ^ + + 0(a , i)2 (3) 


1 


here 
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= a 2 = 2A^ ^ = 1 02 = -l / 4 \ 

m = l, p = 2. K } 

Using equation ( 2 ), the solution of u at any x< = x\ requires the solution of u 
at neighboring values of x<. The solution for a given set of (xi,x 2 ,...x n ) values 
is thus dependent on the solution in some neighborhood of these x values. But 
therefore these neighboring solutions are dependent on their neighbors’ solutions. 
For each Xi this dependency is generally continuous between one or two Xi values 
at which u is definitely known: the auxiliary data. Therefore, the solution process 
starts near values of the x variables for which u is known. One technique, a solution 
by “lines”, is outlined as follows. One of the x variables, say Xi is stepped through 
a set of values. At each value a finite difference equation equal, to some order of 
accuracy, to the original partial differential equation, is set up. The resultant set 
of equations is then solved simultaneously. x 2 may then be varied a small amount, 
Ax 2 , and the process of running through x\ values is repeated to solve another 
set of equations. When x 2 reaches its maximum X3 can then be varied a small 
amount, AX3, and the previous processes repeated. 

The fact that u changes only by changing one x variable makes it convenient 
to view the problem within an n-dimensional coordinate system where each X* is 
represented by an axis. Within this system the set of x values in equation (2) lies 
on a coordinate line 1 of the system. In a two-dimensional system, where x = x\ 

and y = x 2 , u x is calculated from values of x and u along a y coordinate line; u y 

*a coordinate line is described by allowing x, to vary while keeping all other coordinate 
variables constant, a coordinate surface is described by keeping Xi constant while varying the 
other variables. 
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is found from values of y (and u) along an x coordinate line, figure 1 shows the 

J Wi.-wiO Wta** ^ •—.-<* - -e. » .* -- .* vmJ. « uliObii^ta X 0 bii.Om.LX OO 

obvious from the figure that the field is now broken up into a set of points. It is 
at these points that the solution is obtained. 

To facilitate numerical computations it is advantageous to keep Ax and Ay 
constant, for the following reasons. 

• The derivation of the above finite difference equation, (3), is based on a 
constant Ax, allowing for a higher order of accuracy and simplicity of form 
than would otherwise be possible. 

• A constant A leads to the symmetrical array of points in figure 1. These 
points can then be indexed as an array, with points in the x direction being 
indexed by j, and points in the y direction being indexed by k. When one 
boundary value of x , is indexed as a:(l), and the other boundary as x(jmax), 
the following simple relations can then hold 

x(j) = z(l) + ( j - 1) X Ax (5) 

x(jmax) = x(l) + {jmax — 1) x Ax (6) 

Where j = 1,2,. . .jmax. A similar situation exists for y where the index 
runs from one to some maximum, kmax. One can then say that 

X = x(j) = Xj 

y = y(k) = yk ( 7 ) 

u = u(j,k) = U j)k . 

Further, since to each u there associates an x and a y, they can be indexed 
as 

omomt & 
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x(j, k ) = iSi t k j = 1 , 2 ..,. ,jmax 

y(j,k) W yjf ~~ T ‘ -'"***> 

Any point in the field can then be identified by its j,k indices. 


There are three major reasons why the A’s are not usually constant. 

• To reduce errors due to 0(Ax) n and increase resolution of u, Ax can be 
made smaller. 

• To increase efficiency Ax can be made larger in areas where u is changing 
more slowly. 

• Varying y may vary the boundary values of x, so that equation (6) holds 
only if Ax is changed. 

Figure 2 shows the result of varying Ax within the field, in the area labeled 
“increased resolution”. Because of the constancy requirement for a finite dif- 
ference expression, this new point requires additional points along the x and y 
coordinate lines intersecting this new point, u must be solved for at these new 
points, increasing computational time. 

When the boundary does not lie on a coordinate line, as in figure 3, unless the 
boundary is a linear function of y and x, it is not possible to keep Ax’s or Ay’s 
constant in the field. Although interpolation can be employed or different finite 
difference terms used, accuracy tends to decrease, while computational time tends 
to increase. 

One solution is to define a new coordinate system, a computational system, in 
which 


• All boundaries coincide with coordinate lines (or surfaces in 3-d). 
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• A minimum of coordinate lines are used to to describe the boundary. 

• If the new coordinate system’s variables are denoted by £,77 then u = u(£,r}), 
and 

• the finite difference approximations contain terms in £ and tj, and A£ and 
Arj. 

• A£ and Arj, the point spacing on the boundary and interior, are constant. 

• The mapping between coordinate systems must be one-to-one so that 



(9) 

v (z,y), 

(10) 


(11) 


(12) 


• To ensure nonsingularity between corresponding points the Jacobian deter- 
minate (J) must be nonzero. 

• When £,77 are the coordinates of a cartesian system the grid is a rectangular 
mesh of points. 

The requirements lend themselves rather well to using partial differential equa- 
tions as the mapping functions. Existence of 2nd order partials also ensures 
smoothness and continuity. Although it is possible to use conformal mapping 
techniques to map between two-dimensional coordinate systems, conformal map- 
ping does not extend to three-dimensional systems, whereas partial differential 
equations can map between these systems. 
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Figure 1: Portion of boundary points and solution points in a 2-dimensional rect- 
angular coordinate system with Ax and Ay constant. Boundary conforms to a 
coordinate line. 
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Figure 2: Variance in Ax and Ay to increase local solution accuracy and resolution 
with resultant additional number of solution points. 
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Figure 1: Consequent change in Ax and Ay near boundary that doesn’t conform 
to a coordinate line. 
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2 Classification of Partial Differential Equations 


The general form of a second-order partial differential equation in two variables is 

Q'Uxx ~|“ bzixy “I - cttyy d - du x ctiy /w — Q (13) 

where u = u(x,y) and a,b,c,d,e,f, and g axe functions of x and y for a linear 
partial differential equation. If these coefficients contain first partials of u (u x or 
Uy) the equation is termed quasi-linear. 

By an appropriate transformation from x,y to £,77 (13) transforms into one 
of three canonical equations, depending on the value of b 2 — 4®c. 

For b 2 — 4ac > 0 (13) is termed hyperbolic and the canonical equation is 

u ii~ u wi ^ — = 0 (14) 

For b 2 — 4 ac = 0 (13) is termed parabolic and the canonical equation is 

U££ + • • • = 0 (15) 

For b 2 — 4 ac < 0 (13) is termed elliptic and the canonical equation is 

Utf + u VT] + • • • = 0 (16) 

Where the • • • refers to terms in x, y, tt, U(, and u v . Each type of equation defines 
what auxiliary data in u, or one of its derivatives, is required for the problem to 
be “well-posed” for solution. This also directly influences what solution procedure 
may be used. If both x and y are spatial functions the auxiliary data is a boundary 
condition. If x or y represents time the auxiliary data represents initial conditions, 
and possibly final conditions. Both hyperbolic and parabolic equations march out 
from the auxiliary data to an unspecified boundary (or final ) condition of u. 
For hyperbolic equations, however, the solution at a point within the domain of 
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integration is a function only of a certain portion of the given data. For parabolic 
equations, she solution at a point is dependent on the entire auxiliary uata. ise 
elliptic equation requires that the auxiliary data in u surround the domain of 
integration. That is, u for the entire boundary of the domain of integration must 
be given. The solution is then found for u within the domain. 
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3 Transformation Metrics 


Given a general curvilinear coordinate system denoted by coordinates £ and 77 , 
with a one-to-one point mapping between it and a cartesian system denoted by x 
and y, that is 

£ = ( 17 ) 

77 = T)(x,y) (18) 

as well as 


« = ( 19 ) 

y = ( 20 ) 

then by the chain rule, 

d£ = £* dx + tydy (21) 

dx = + x v drj (22) 

dy = ytdt + yndrj. (23) 


Substituting the values for dx and dy of equations (22) and (23) into equation 
(21) produces 

d£ - tx(xid£ + x v dr)) + £ y (^d£ + y v dr)). (24) 


Rearanging, 
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Since £ is not a function of tj the terms within the second set of parentheses 
must equal zero and so the terms in the first set of parnetheses must equal 1. The 
resulting set of equations 


+ £y2/£ — 1 

£x x t) "t“ ^ylfri = 0 


produces 


tx 

€y 


yy_ 

\j\ 

-X, 

\J\ 


( 26 ) 

( 27 ) 


where 


\J\ = HVn ~ x v Vi 


is called the Jacobian determinant. The same procedure can be applied to r/ and 
the results are 


ORfct&U. i'£ 

of. poor quality 

4 Elliptic Grid Generation 

The simplest, and most commonly used partial differential equation for grid gen- 
eration is Laplace’s equation [3], [2] 


V It — “J“ Uyy — 0 


(30) 


a homogeneous elliptic equation. 

Replacing u in turn by the transformation coordinate variables £ and 7 / pro- 
duces a system of equations 


where 



£xx + €yy 
Vxx + Vyy 



(31) 



When £,7/ are considered as the independent variables, equation (31) can be 
written, using the transformation metrics, as 


aIR a - 2 (3IR ir] + 7 IR vr , = 0, (32) 


where 


R 

a 

P 

7 

I 



2 , 2 


+ yty v , 


x\ + 2/f , 
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Using central difference forms for the first and second derivatives (s< 
pendix), noting that A£ = A 77 = 1 , and using indexing notation, one 
difference form of equation (32) is 


qti 

Hk (2A(y 


2b, a 


S^R 


(2Af)(2A 77) 


+ c i*I\ 


(2A v y 


= 0 , 


where 


a j,k 

bj,h 

c j,k 


/ ^t] X j,k x2 , / ^yVj,k \ 2 

K 2A 7] } ~ r{ 2Arj ) 

/ ^Vj,k \ / ^yVj,k X 

V 2A£ A 2At/ ' 1 2A£ A 2At; ' 

( ^ X j,k \2 , / ^Vj,k \2 

V 2A£ ; ri 2Af 


After clearing the denominators and expanding, equation (33) becomes 


a j,kI(Rj+i,k 2Rj,k 4 " Rj— I,*) 

- 26 J) fci’(^ + i i fc_i 4 - i,*.— i) 

2,Rj t k 4* Rj,k—i) = 0. 


Gathering all il terms in k to the left side, this now becomes 


Aj,kI(Rj-i,k) 4- Bj t kI(Rj,k) 4- C jtk I(R j+1) k) — Gj fk 


e *ap- 
finite 

(33) 

(34) 

(35) 

(36) 


(37) 

(38) 
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where 

Aj,k — C jtk = aj t k 

Bj,k — — 4- c j,fc) 

Gj,k = %bjjkI\Rj+ i,fc+i — -Rj— i,fc+i “ -Rj+ijt— i 4“ -Rj— i,A-i ) 

~ c j>kI{Rj,k+ 1 H - Rj,k—i) 

This system can be solved by the “line” solution technique mentioned earlier: 
by keeping k constant (starting at k=l) an equation for each Rj t \ (from j=l to 
j=jmax) is generated. The resultant set of equations is then solved simultaneously 
for R. The index variable k is then incremented and the process, of generating a 
set of equations for all j, is repeated. The overall process of incrementing k and 
solving a set of equations continues until a final k=kmax value is reached. Because 
the values of all x and y are initially unknown, except at boundary points, initial 
guesses for them must be made. These initial x and y values are then used to 
generate the values for A,B, and G. As new values of Rj t k ( %j,k and yj^ ) are 
found, they are used to solve the succeeding set of equations for the next set of R 
values. 

equation (38) is a system of 2 equations, one for each unknown coordinate 
variable x, and y. Substituting x for each instance of the vector R the equation 
becomes 


Aj, k Xj-i,k + Bj,kXj,k + C jt kXj + i t k = 

2bj,k( x j+l,k + 1 ~ *i-l,*+l “ ®i+l,Jb-l+j-l,k-l) 

- Cj,k(j,k+l+j,k-l)- ( 39 ) 


The complete set of equations, for constant k, from j=l to jmax, can be put 
in matrix form as 
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( 40 ) 


( B2 ,u Cz,k 0 . . . 0 

Ai t k . Bzj, C* ,k 0 ... 0 

0... ••• : 

0 ...0 Aj-i t k Bj t k Cj + i'k 0...0 

• • • 

• • • 

• • « 

0 • • • 0 Aj max —2,k Bj max —i'k 


/ *2 ,fc ^ 


1 Gx 2 ,fc ) 

* 3 ,fc 


Gx 3l k 

X i.k 

~ 

Gxj'k 

\^jmax—l,kj 


\Gxj max _i i f e J 


where each Gx^k equals the right-hand-side of equation (39). A similar set of 
equations are solved for y at each k level. 

Having run through all j and k values an improved set of x,y values is obtained. 
This improved set can then be used to obtain an even better set of values by 
repeating the procedure just described. This iterative process stops when at some 
iteration the change in x and y from a previous iteration is smaller than some 
given amount. 

The number of iterations can be reduced significantly by using a technique 
known as successive line over-relaxation, [6] which is essentially a forecasting tech- 
nique. The latest values of x and y are compared with their previous values and 
some fraction of their difference is added to the latest value, to predict a better 
value for the x and y’ s. For example, after solving the set of equations along a 
line using Thomas’ algorithm, and obtaining new values for x(j,k) at each k level, 
where n refers to the iteration number, the following modification can be applied. 


x lt = x lt' + u ( x lk - x lt') 

where n references the last iteration of x and n-1 references the previous value. 
u> controls the amount of change made. When u> lies between one and two the 
technique is known as over-relaxation. When u> is between zero and one the 
technique is known as under- relaxation. Under-relaxation is useful when values 
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of the dependent variable oscillate for a given set of values of the independent 
variables, that is, oscillation with respect to iteration 

The main advantage of Laplace’s equation in grid generation is that it serves 
to smooth grid lines. This is because it contains within it the maximum principle: 
there are no maximums or minimums of u within the domain of integration. That 
is, u attains its highest or lowest values on the boundary of the domain, u(x°, y ° ) is 
an average of u in the neighborhood of x° and y°. This serves to reduce disparities 
in u for neighboring points. This can be shown by noting first that, within the 
domain of integration, the dependent variable cannot have a maximum or mini- 
mum. For there to be an extremum the first partial derivatives of the dependent 
variable must each go to zero; the second partials must have the same sign and 
each must be nonzero. This would mean u xx + u y y ^ 0 which is not the case here. 
Since u does not have an extremum within the boundary it must be bracketed by 
values of u on the boundary. 

Figure 4 is an initial guess for a grid about a body. This grid was generated 
by placing points evenly on a circular outer boundary and joining them to points 
on the inner boundary, the body. Then points were evenly spaced along these 
connecting lines to generate the nodes. Figure 5 shows the results of 100 iterations 
of the elliptic solver. There has been a general smoothing of lines throughout. 
Figure 6 is a detail of the initial grid in figure 4 showing grid breakdown near one 
of the convex points. Figure 7 is this same area after 100 iterations of the elliptic 
generator 

The main disadvantage to using this system of equations is that there is no 
control over grid point spacing. For each given set of boundary values there is 
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only one final distribution of interior points. An alternative is to use the nonho- 
mogeneous form of the elliptic equation, known as Poisson’s equation 


V 2 £ = P (41) 

V 2 7/ = q (42) 

In such equations p and q can be used to modify the strong tendency of the 
Laplacian to equalize all cell volumes [7]. Another alternative takes the fully 
converged grid and re-spaces the points along £ lines in order to cluster points 
near the body. [4] 
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Initial Grid 



Figure 4: tau=.5, cam=0.8, Initial Grid 
File: s5t8.00ps 
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The Effect of the Elliptic Solver 



Figure 5: tau=.5, cam=0.8, pass 100 
file: s5t8.10ps 
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Initial Grid 


i 



Figure 6: tau=.5, cam=0.8, detail of initial grid of figure 4 
file: s5t8.00ptps 
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The Effect of the Elliptic Solver 



Figure 7: tau=.5, cam=0.8, pass 100, detail of final grid of figure 5 
file: s5t8.10ptps 


22 




5 Hyperbolic Grid Generation 


Steger [8] developed a grid generation scheme based on a set of equations which 
together form a hyperbolic system. These equations are expressions of constraints 
made on the grid generating system. These constraints are as follows: 

• Orthogonality is explicitly specified, 

V£ • V?7 = £ x 7) x + £yT)y = 0. (43) 

• Grid cell volume (i.e. area in two dimensions), V, is user defined. Since, in 
the numerical implementation A£ = A rj = 1, 

dx dy = (x^yrj - x v y^)d^ dr] = J~ l d £ drj (44) 

d£ dr] — ^2/ — J dx dy (45) 

Thus the Jacobian should be nonsingular and approximately equal to the 
physical cell volume, Ax Ay. 

The resultant system is thus represented by 

CxVx + £yVy = 0 (46) 

CxVy ~ £yVx — J ( 47 ) 

(48) 

or, after a change of variables, again using the transformation metrics, 

**«i» + 3fc Vy = 0 

x iVrt ~ x vV( = J~' = V 
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(49) 

(50) 



Equations (49) and (50) comprise a system of nonlinear partial differential 
equations with initial data in tj. Through a process known as local linearization, 
shown below, it is found that this system is hyperbolic and thus can be marched 
in r] from initial data along the body, rj = 0. 

If the substitutions x = x° + x and y = y° + y are made, where x° and y° 
represent a known nearby state, and x and y are of the order of Ax and Ay, 
local linearization can be performed by making these substitutions for x and y in 
equations (49) and (50). Starting with the first equation, 


X(X V = (z° + x)t)(x° + x) v (51) 

= x°^x° v + X^Xr, + X° V X£ + X£X V (52) 

= xl^n + x^x v + x\x^ (53) 

= X^Xrf + xl(x£ - X \x Q r (54) 

Similarly, the remaining terms in the two equations become 

vm = ylvn + ylvi-ylyl ( 55 ) 

XiVr, = yftt + XtVr, - x\yl (56) 

x v y( = y^x v + x° v yt - x^yl (57) 


Using these substitutions equation (49) becomes 

+ A x i - x l x v) + (y^vv + vlvt - y\yl) = 0 ( 58 ) 

But since X(X V + y ^ = 0 this reduces to 

x 0 ^ + z°xt + y° c y v + yfy = 0 ( 59 ) 
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or, in matrix notation 



Doing the same for equation (50) 

(vl x t + x \yr, - x \vl) - (y^v + x ° v yt - = v (6i) 

Since by equation (50), — X ^V^) = V° the equation becomes 

yl x t + x hv -y&v- x lyt = v + v ° ( 62 ) 


or 



Combining equations (60) and (63) 




( 64 ) 


which can be represented as 


AR^ + BRrj = /• 


( 65 ) 
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If B, whose determinant is the Jacobian transformation, is non-zero, all cell 
volumes are non-zero, and its inverse, B _1 , exists. Equation (64) can be then be 
transformed into 


B-'AR ( + £, = B~ l f. 


( 66 ) 


It has been shown [1] that this type of equation is hyperbolic and marches in 
the rj direction if B~ X A has real, distinct eigenvalues, or if B~ x A is a symmetric 
matrix. Expanding B~ x A and performing the indicated multiplication produces 


B~ l A = 


( 4 ~yt ) 

U? 4 J ( 4 y°r, ) 
44 - l yl -4 ) 

( 44 - yK 4$ + 44 \ 
[44+44 44-44J 
(4) 2 + (yl) 2 


( 67 ) 

( 68 ) 


which is symmetric, fulfilling the latter requirement. If A represents the eigen- 
values the eigenvalue equation is 


( - y°y°) x) 2 

( ylx° v + x\yl \ 

l K) 2 + (y?) 2 ) 

v(*?) 2 + (y«) 2 / 


( 69 ) 


or 


> Yl ± ( x K + gM) 

* K ) 2 + ( y ?) 2 

Thus, since V°, the inverse Jacobian, is non-zero the eigenvalues are real, and 
distinct when is non-zero 
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5.1 Numerical technique 

The numerical solution of equation (66) can be carried out using the line inversion 
mentioned earlier in the introduction and in the section on the elliptic generator. 
The main difference here is that only one pass is performed in generating the 
grid. Unlike the elliptic generating system, which improved an initial grid, and 
therefore was dependent on the quality of that grid for convergence to a solution, 
the hyperbolic system generates a grid from the known boundary conditions. 

Replacing the partial derivatives in equation (65) by finite difference operators 
produces 


-1 a I Vtfi _ 


B~ A 


= B~'f 


2A£ At) 


( 71 ) 


where all terms are second order accurate except for the first order replacement 
of R,,. 

Since, in the numerical implementation, A£ = A 77 = 1, the expansion of 
equation (71) is 


B~ X A 


{Rj+l,k+l Rj—l,k+l) + Rj,k + 1 Rj,k — B f (72) 


Rearranging to gather terms involving k + 1 on one side produces 


^j,kRj—l,k+l “1“ ^j,kRj,k+lCj } kRj+l,k+l — RH S 


( 73 ) 


where 


Hk 


B~ l A 

2 
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( 74 ) 



H 

II 

(75) 

B~ l A 

(76) 

c ),t - 2 

RHS = B-'f + R hk 

(77) 


The coefficients are all evaluated at the same index values of j and k, where 
x and y are known. The partials with respect to £ are straightforward central 
differences. The rj partials can then be obtained by solving for them in equations 
(49) and (50) 


Vv 


yj V 

(^) 2 + M 2 

XjV 

M 2 + fo ) 2 


(78) 

(79) 


Equation (73) generates the set of equations in j, from j=l to j=jmax. These 
equations are solved, as before, using Thomas’ algorithm. 

The results of using this algorithm, designated hyg2d, can be seen in figures 8 
through 11 for bodies of increasing concavity 2 . Note that as the bodies increase 
in concavity the radial grid lines start to cluster, and finally, the grid breaks 
down, figure 12 and 13 show the marching out of the grid just before and at this 
breakdown. Figure 14 is a detail view of one point of the body. 


2 See the section Notations and Definitions and figures 69-72 for more information on concavity 
and how it is meant here. 
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The Effect of Camber (concavity) on the Hyperbolic Solver 



Figure 8: tau=.5, cam= 0 
file: h50y0.ps 
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The Effect of Camber (concavity) on the Hyperbolic Solver 



Figure 9: tau=.5, cam= 0.40 
file: h50y40.ps 
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The Effect of Camber (concavity) on the Hyperbolic Solver 



Figure 10: tau=.5, cam= 0.80 
file: h50y80.ps 
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The Effect of Camber (concavity) on the Hyperbolic Solver 



Figure 11 : tau=.5, cam= 1.0 
file: h5yl0.ps 
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Marching Out of Hyperbolic Generator 



Figure 12: tau=.5, cam=1.0, just before cavity-caused breakdown 
File: h5yl0.kl5ps 
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Marching Out of Hyperbolic Generator 



Figure 13: tau=.5, cam=1.0, breakdown due to strong concavity 
file: h5yl0.kl6ps 
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The Effect of Camber (concavity) on the Hyperbolic Solver 



Figure 14: tau=.5, cam= 1.0 detail 
file: h5y!0.ptps 
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6 Combining Elliptic and Hyperbolic Equations 

Each of the grid generation schemes discussed, the elliptic and the hyperbolic, 
has a profile of advantages and disadvantages. The elliptic, through its smoothing 
action is able to generate a grid about more complex bodies than the hyperbolic. 
But it requires an initial grid to start from; it must iterate to convergence for 
a solution; Unless source terms are added, there is little control over point dis- 
tribution; these source terms are not easy to use. The hyperbolic, on the other 
hand, creates a completed grid in the order of time it takes the elliptic to make 
one pass. This particular hyperbolic scheme has, additionally, orthogonality con- 
straints, and user-control of grid cell volumes, providing a strong control over 
point distribution. Its chief disadvantage is that unlike the elliptic it breaks down 
sooner in areas of strong concavities or convexities. 

A scheme that uses the best of both would be better than either. Steger [5] 
has developed such a grid generation scheme. His procedure is to add a scalar 
multiple of the elliptic system, equation (32), to the right side of the hyperbolic 
system, equation (66), producing 

HYPERBOLIC : 

B-'ARt + Rr, = B~ l f ( 66 ) 

fl X ELLIPTIC : 

0 = - ‘ZpRtn + 7 Rrm) ( 32 ) 

COMBINED : 

B^ARt + R v = B~ l f + fi(aRtt - 2 f3R ir) + 7 R^) ( 80 ) 

The control parameter fi can then be used to emphasize the hyperbolic or el- 
liptic system. For small /z the hyperbolic dominates and for large /z the elliptic 
dominates. 
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Viewed as a hyperbolic system the constraint of orthogonality is now relaxed 
by the strength of the elliptic terms when they sum to nonzero. Loosening this 
constraint helps avoid grid breakdown. In the volume equation the elliptic terms 
serve to smooth differences between cell volumes. Viewed as an elliptic system, the 
Laplace equation has been transformed into a Poisson equation with the hyper- 
bolic terms serving as the source terms. These source terms allow concentration 
of £, and r) lines. Which side dominates is controlled by ft. When it equals zero 
the system is purely hyperbolic. As it increases from zero the right side of the 
equations become more dominant until fi is over about 3,000, when the system is 
effectively elliptic. This parameter can be made a function of spatial position, in 
order to, for instance cluster lines near the body. It can also be made a function 
of gradients in the flow field so that, in an area of high gradients, fi can be set 
low, turning on the hyperbolic’s volume control, to increase point density in this 
area. In areas of grid breakdown fi can be set high, turning on the elliptic, to 
repair the breakdown. 

6.1 Solution Procedure 

Replacing the partial derivatives in equation (6) by finite difference operators 
produces 


B-'f + V [«#* 


(A£) : 


-2 (3 


S ( S V R 


. & 

+ Tr^ 


(2AO(2Ai,) '(A ,)» 


] (81) 


where all terms are second order accurate except for the first order replacement 
Of Rr,. 
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Since, in the numerical implementation, = A-rj = 1, the expansion of 


equation (81) is 


B~ l A - - 

— — ( R j+l,k - Rj-l,k) + Rj,k - Rj,k - 1 = 


B V 4- va(R j+lt k - 2R jt k + Rj-i,k) 

+ ^-{Rj+imi ~ Rj-ijt+i ~ Rj+i,k-i + -Rj+i,fc-i) 

+ ^(Rj'k+i — 2Rj,k "t - Rj,k- i) 


(82) 


Rearranging to gather terms involving k on one side produces 


Q'j,kRj—i i k “ 1 “ bj,kRj y k Cj-\-i,kRj,k — RH S 


(83) 


where 


Hk 

b 3,k 

c j,k 

RHS 


B~ l A 
2 

1 + 2fi(a + 7 ) 

B-'A 

— 

Bj,k—1 “I" 7 fl{Bj,k+l 4" Bj k—l ) 


(84) 

(85) 

(86) 


— 9 (-fy+i,*+i ~ Rj-iM 1 — Rj+ijt-i + -Rj+i,t-i)(87) 

( 88 ) 


and all values in equations (84) - (88) not specifically indexed have index values 
of j and k, that is, they are evaluated at the points x = Xj t k, y = yj,k • 

Equation (83) can now be solved similar to the elliptic equations. Starting with 
k=l, letting the j index run from 1 to some jmax, and keeping the k index constant, 
a set of equations are generated. The left-hand coefficients form a tridiagonal 
matrix, each element of which is a two-by-two array. The set of equations is 
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solved using a modified form of Thomas’ algorithm. The index variable k runs 
from 1 to some number kmax and therefore kmax sets of equations must be solved. 
Once these are solved the whole process is reiterated using the results of the last 
iteration, the resultant grid, as the source for generating new values for the right- 
hand side and the coefficients, of equation (83). The iteration process halts when 
grid changes are below some threshold value. 

6.2 Testing Procedure 

Since the elliptic can repair breakdowns in the grid (due to the hyperbolic), and 
since breakdown is a function of the number and strength of convexities and 
concavities, it was decided to test this new solver by using bodies of increasing 
concavity and/or convexity . 

Although the hyperbolic generator breaks down in areas of either great con- 
cavity or great convexity, breakdowns due to the latter tend to be more localized 
and simpler to deal with. Uneven point spacing along the body on either side of 
the apex of the convexity causes the grid line radiating from the apex to swing 
over. This causes circumferential lines to cut through the body. But the number 
of nodes involved are generally limited to a few in either the circumferential (£) or 
radial ( 77 ) direction and can be handled by the fortified approach of Van Dalsem 
[9]. 

The breakdowns caused by a concavity tend to be more global as the strength 
of the concavity, and the percent body surface the concavity takes up, increase. 
Therefore, it was decided to approach the problem by examining bodies containing 
concavities of increasing strength. 
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6.3 The Unmodified Algorithm 


Implementation of the elliptic-like algorithm just described requires an initial grid. 
The main requirement placed on this grid is that its outer boundary not self- 
intersect or intersect with interior grid lines. The hyperbolic generator was used 
to generate the initial grid of figure 15. This generator offers the advantage of being 
closely related to the hyperbolic part of the combined generator. The similarity 
can be used to avoid problems such as “collapsing” grid lines, if lines that tend to 
fall back towards the body rather than away from it with succeeding iterations. 
This defect can occur when the initial grid has different spacing in the 77 direction 
than the hyperbolic component provides. Here p is set low near the body and 
increases geometrically in the radial direction. Thus the hyperbolic is enforced 
close to the body allowing clustering of points, while the elliptic is more in force 
deeper into the field out to the outer boundary, smoothing the sharper curvature 
of grid lines of the initial grids. Figure 16 shows the grid after 100 iterations. 
Note that because the outer boundary the clustering on it cannot be dealt with. 

An alternate initial grid, having a strong concavity and the final grid drawn 
from it are shown in figures 17 and 18. Figure 19 shows detail of the latter 
grid, at a convex point. This initial grid was generated by specifying a circular 
outer boundary and spacing the points on it evenly. Lines are drawn connecting 
these points with points on the body. Points are then evenly spaced along these 
connecting lines to form the grid. Here the outer boundary point spacing was 
fixed to avoid the clustering of the hyperbolic system. For high concavities it 
seems to generate a reasonable grid. Because the outer boundary point spacing 
is not a strong function of the inner boundary, their connecting lines (£ lines) 
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sometimes have to make large sweeps. As shown in figure 19 orthogonality can 
break down. Other drawbacks are that there is no way of knowing which points 
are best connected to each other and, eventually a body shape is reached at which 
this algorithm breaks down. 

6.4 Floating the Boundary 

The major problem with the original algorithm is that it found a solution in 
an elliptic way. That is, the boundaries were kept fixed while the inner field 
points were adjusted. The first modification to the algorithm above, and used in 
subsequent modifications, was to float the boundary. This was achieved by turning 
the hyperbolic on at the outer boundary to generate a new outer boundary. This 
technique is more amenable in areas where there is no physical outer boundary, 
such as for airfoils, areas that use hyperbolic generators to form grids. 

Once the outer boundary is floated initial grids containing defects in their 
outer boundary, such as those generated by hyg 2 d, can be used. Such defects 
include boundaries that intersect other parts of the grid. The smoothing action 
of the elliptic, coupled with the boundary generating action of the hyperbolic can 
then repair defects in the initial grid. Figures 20 to 33 are a series of grids showing 
the effect of floating the boundary. 

The drawback to this latest technique is that at some point a body is reached 
whose concavity makes it too difficult to generate a grid about. The generator 
then fails to improve the grid. Even before this point is reached, though, the range 
of values \i can take on narrows, making it difficult to choose a set of values for it. 
Here fi is just a function of 77 . For bodies with weak concavities it can take on any 
value at any 77 . But as concavity grows grid improvement becomes very sensitive 
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to the value of /x, to the point that it becomes very difficult to find values for /x 
that will enable a solution to converge. 

6.5 Using Compressed- Cell Initial Grids 

The problem with the previous technique was in the initial grids used. They 
contain too many defects and slow the convergence procedure. If cell volumes are 
extremely small the distance from inner to outer boundary makes intersections 
of grid lines more difficult. This led to modifying the algorithm to handle a 
different kind of initial grid, one with fewer defects. An initial grid of compacted 
cell volumes was used. During each iteration all cell volumes were then increased 
a fixed percent of their total volume. This had the benefit of starting with an 
initial grid which had very few, if any, defects and smoothing out defects as they 
developed, rather than starting with a grid containing many defects. Defects were 
generally isolated then, at or near the outer boundary. Each iteration was both 
a growth of the grid, somewhat like the hyperbolic, and also each iteration was a 
smoothing out in the interior, by the elliptic. Since the outer boundary is purely a 
product of the next lowest level, which was constantly being improved, this made 
for a very robust generator. Compared to the previous method this technique 
allows for a wider range of values for /x given the same body as before. Figures 
34 to 45 are examples of this technique for bodies of high concavity. Cell volumes 
have been set to grow about fifteen percent per iteration. Figure 34 contains the 
initial grid but, since its cell volumes are compressed the cells can only be seen 
at at several magnifications (figure 35). Figures 36 to 40 show the grid expanding 
during successive passes of the combined solver. Figure 41 is a close-up of the the 
grid at the 200th pass. Figures 42 to 45 shows that the technique can be applied 
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OF POOR QUALITY 

to even deeper cavities. 

6.6 Curtailing Outer Boundary Movement within a Cav- 
ity 

There are two main drawbacks to using the technique just discussed as it stands. 

• A point is reached in which the cavity is too deep for this technique to work, 
see figure 46 

# Although a grid with no breakdown is possible for strong concavities, there 
tends to be less control over point distribution within these areas. 

The technique in the previous section uses just the hyperbolic component to 
generate a new outer boundary with each iteration. The hyperbolic generates 
the points on this boundary, the rjkmax line from data on the next lower section 
of the grid, the T)kmax-i line. By the orthogonality constraint these points are 
normal to the Tfkmax-i line. As the cell volumes increase this causes boundary 
self-intersection within a cavity with some depth. Because the interior is derived 
primarily elliptically and therefore a function of both boundaries there is a ten- 
dency therefore for grid breakdown near the outer boundary. But, since the inner 
boundary is okay, there is, by the maximum principle less breakdown of these 
interior grid lines than the boundary line. At each iteration these areas of break- 
down are improved by the elliptic through communication with inner grid points. 
These improved areas are then used to generate the latest outer boundary. The 
fact that the numerical technique starts at the inner boundary may be part of the 
reason why the elliptic smoothes these areas so well. 


43 


The result of the above tendencies is that within the cavity the outer boundary 
moves, at each iteration, largely in the 77 direction and out of the cavity. Thus in 
strong concavities the ratio of cell side lengths can be very high. See figure? And, 
while volume control normally gives control over point distribution, in this case it 
is overidden by the smoothing action of the elliptic coupled with the hyperbolic’s 
regeneration of the outer boundary, forcing the boundary out of the cavity. Within 
the cavity accuracy and resolution are thus lower in the 77 direction because of the 
resultant spacing. It was thotight^aitnfksk that' point -density was lower, in vthe 
cavity, contributing to the increased cell aspect ratio, but analysis confirmed that 
for a given range of 77 point density is about the same within and without the 
cavity. This conforms with the fact that cell volumes had been functions of 77 
only. 

In and near an area of a concavity the gradients of the flow variables can be 
expected to be higher than average. Therefore A^ and should actually be 
smaller than elsewhere to resolve the flow and provide a stable solution. There 
then appears to be a conflict between the need to maintain a small cell size and 
at the same time the need to cover the flow field out to several body diameters. 
It was therefore decided to modify how the outer boundary is generated within 
the cavity. A point on the boundary is chosen to serve as a hinge point. As the 
boundary points move away normal to the i]kmax-i line their movement is modified 
to some specified degree by a forcing function that brings corresponding points, 
on each side of the hinge point, together. When a set of points meet they are no 
longer allowed to float. Their corresponding cell volumes are then held constant 
through later iterations. Figures 47 to 62 show the results of using this technique 
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on the grid that caused breakdown of the previous technique. There is The result 
is that within the cavity the boundary overlaps itself along a line of symmetry 
of the cavity and it is now possible to generate a grid for deeper cavities also a 
higher point density within the grid. 





Using the Combined Solver with a Fixed Outer Boundary 



Figure 15: tau=.5, cam= 0.8, initial grid 
file: eh5y8.00ps 
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Using the Combined Solver with a Fixed Outer Boundary 



Figure 16: tau=.5, cam= 0.8, pass number 100 
file: eh5y8.10ps 
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The Effect of Camber (Concavity) on the Combined 
Elliptic-Hyperbolic Grid Generator 



Figure 17: tau=.5, cam=2.0. Initial Grid 
file : s5t20.p0p.stmp 
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The Effect of Camber (Concavity) on the Combined 
Elliptic-Hyperbolic Grid Generator 



Figure 18: tau=.5, cam=2.0, pass 100 
file: seh5t20.10ps 
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Figure 19: detail of figure 18 
file :seh5t20.10ptps 
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Grid Improvement by Floating the outer Boundarv 



Figure 20: tau=.5, cam=1.0, initial grid 
file: ehf510.00ps 
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Grid Improvement by Floating the outer Boundary 



Figure 21: tau=.5, cam=1.0, pass number 1 
file: ehf510.01ps 
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Grid Improvement by Floating the outer Boundary 


Figure 22: tau=.5, cam=1.0, pass number 2 
file: ehf510.02ps 
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Grid Improvement by Floating the outer Boundary 



Figure 23: tau=.5, cam=1.0, pass number 3 
file: ehf510.03ps 
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Grid Improvement by Floating the outer Boundary 



Figure 24: tau=.5, cam=1.0, pass number 4 
file: ehf510.04ps 
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Grid Improvement by Floating the outer Boundary 



Figure 25: tau=.5, cam=1.0, pass number 5 
file: ehf510.05ps 
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Grid Improvement bv Floating the outer Boundai 



Figure 26: tau=.5, cam=1.0, pass number 6 
file: ehf510.06ps 
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Grid Improvement by Floating the outer Boundary 



Figure 27: tau-.5, cam=1.0, pass number 7 
file: ehf510.07ps 
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Grid Improvement by Floating the outer Boundarv 



Figure 28: tau=.5, cam=1.0, pass number 8 
file: ehf510.08ps 
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Grid Improvement by Floating the outer Boundary 



Figure 29: tau=.5, cam=1.0, pass number 9 
file: ehf510.09ps 
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Grid Improvement by Floating the outer Boundary 



Figure 30: tau=.5, cam=1.0, pass number 10 
file: ehf510.1ps 
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Grid Improvement by Floating the outer Boundary 



Figure 31: tau=.5, cam=1.0, pass number 20 
file: ehf510.2ps 
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Grid Improvement by Floating the outer Boundary 



Figure 32: tau=.5, cam=1.0, pass number 50 
File: ehf510.5ps 
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Grid Improvement by Floating the outer Boundary 



Figure 33: tau=.5, cam=1.0, pass number 100 
file: ehf5 10.1 Ops 
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Figure 34: tau=.5, cam=2.0, full figure 

file: hg5y20.ps, initial dseta=le-5, final dseta=le-3 
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TnitifU r '-~ A Compressed Cells 



Figure 35: tau=.5, cam=2.0, detail of initial grid near convex point 
file: hg5y20.ptps initial dseta=le-5 final dseta=le-3 
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or poor quality 


Cell Volume Gtwth Counted with Ffnatiwt Boundarv 
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Figure 37: tau=.5, cam=2.0, 50th pass 
file: hg520.5elxps 
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Cell Volume Growth Coupled with Flontm^ Boundary 



Figure 38: tau=.5, cam=2.0, 80th pass 
file: fg520.8elxps 
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Cell Volume Growth Coimled with Floating Bourdarv 



Figure 39: tau=.5, cam=2.0, 100th pass 
file: fg520.10elxps 


70 





Cell Volume Growth Coupled with Floating Boundary 



Figure 40: tau=.5, cam=2.0, 200th pass 
file: fg520.20elxps 
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Cell Volume Growth Coupled with Floating Boundarv 



Figure 41: tau=.5, cam=2.0, 200th pass, detail near convex area 
file: fg520.20elxpts 
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Initial Grid with Co m pre ss ed Cells 



Figure 42: tau=.5, cam=4.0. Full Figure 

file: hg5y40.ps initial dseta=le-5 final dseta=le-3 
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Cell Voka^Gn^thiCouplediWitb^oatrag Boundary 



Figure 43: tau=.5, cam=4.0, 50th pass 
file: fg540.5ps 
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Cell Volume Growth Coupled with Floating Boundary 



Figure 44: tau=.5, cam=4.0, 100th pass 
file: fg540.10ps 
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Cell Volume Growth Coupled with Floating Boundary 



Figure 45: tau=.5, cam=4.0, 200th pass 
file: fg540.20ps 
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Cell Volume Growth Coupled with Floating Boundary 
Failure of Technique Due to Depth of Cavity 



Figure 46: tau=.5, cam=6.0, 500th pass 
file: fg560.50ps 
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Of POOR QUALITY 


Forming a Cut lira Cavity 
Initial Grid 



Figure 47: tau=.5, cam=6.0, initial grid 
File: fgc560.00ps 
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Forming a Cut In a Cavity 



Figure 48: tau=.5, cam=6.0, initial grid, detail of cavity 
File: fgc560.00ptps 
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Fo r mi n g a Cut In a Cavity 



Figure 49: tau=.5, cam=6.0, pass 20 
file: fgc560.2ps 
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Forming a Cut In a Cavity 



Figure 50: tau=.5, cam=6.0, pass 20 , detail of cavity 
file: fgc560.2ptps 
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Forming ?, Cut In « Cavity 



Figure 51: tau=.5, cam=6.0, pass 30 
file: fgc560.3ps 
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Forming a Cut In a Cavity 



Figure 52: tau=.5, cam=6.0, pass 30 , detail of cavity 
File: fgc560.3ptps 
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Forming a Cut In a Cavity 



Figure 53: tau=.5, cam=6.0, pass 40 
file: fgc560.4ps 
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Forming a Cut In a Cavity 



Figure 54: tau=.5, cam=6.0, pass 40 , detail of cavity 
file: fgc560.4ptps 
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Forming a Cut In a Cavity 



Figure 55: tau=.5, cam=6.0, pass 50 
File: fgc560.5ps 
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Forming a Cut In a Cavity 



Figure 56: tau=.5, cam=6.0, pass 50 , detail of cavity 
file: fgc560.5ptps 
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Forming a Cut In a Cavity 



Figure 57: tau=.5, cam=6.0, pass 80 
file: fgc560.8ps 
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Figure 58: tau=.5, cam=6.0, pass 80 , detail of cavity 
file: fgc560.8ptps 
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Forming a Cut In a Cavity 



Figure 59: tau=.5, cam=6.0, pass 100 
file: fgc560.10ps 
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Forming a Cut In a Cavity 



Figure 60: tau=.5, cam=6.0, pass 100, detail of cavity 
file: fgc560.10ptps 
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Forming, a Cut In a Cavity 



Figure 61: tau=.5, cam=6.0, pass 300 
file: fgc560.30ps 
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Figure 62: tau=.5, cam=6.0, pass 300 , detail of cavity 
file: fgc560.30ptps 
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7 Grid Breakdown 

rormm? n t lut err n t.'avttv 

Since it is required that the mapping between real space and computational space 
be one-to-one, that is 


R = R{C) 

(89) 

♦ft* 

T o 

II 

T o 

(90) 


The intersection of grid line segments, other than at nodes, would violate this 
requirement for the following reason. The intersection of these grid lines is an 
intersection of lines connecting points which have adjoining index values. For 
instance one line segment may join P(j,k) with P(j+l,k) the other may join P(i,l) 
with P(i,l+1). Since £ and tj are linear functions of these indices and the indices 
are different in each set, the range of values of £ and tj are different between one 
set of points and the other set. Thus at the point of intersetion in the r space 
coordinate system two values of G exist for the one value of R, violating the 
one-to-one requirement. 

Since the only grid line that can be assumed satisfactory is a given boundary 
line, a check for line intersection should most logically start there and work towards 
the other boundary. But a more elementary check should be to first check whether 
the grid points nearest the boundary, and therefore the lines connecting them to 
the body, are outside the body, as they should be. If k indexes points in the radial 
or tj direction, and j indexes points in the circumferential or ( direction, then these 
first interior grid points can be noted as P(j,k) = P(j,2) = Pj 2 . Pj , 2 joins the 
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body at the point Pj \ which is between the points Pj± i,i* 



Figure 63: Portion of boundary and first row of grid cells next to boundary. Here 
k=2 

If one defines following vectors 


II 

to 

1 

t— L 

(91) 

n = fi- 1.1 - r jt 1 

(92) 

= Tj+ 1,2 - r jA 

(93) 

then, for j increasing clockwise around the the body and k increasing radially 
away, it can be shown that if either of the following sets of conditions are met 

for fj X r 2 > 0 

(94a) 

r 1 x r t > 0 

(946) 

and 


ft XT 2 > 0 

(94c) 

or 


forri x r *2 < 0 

(95a) 

either 


f i x f t t > 0 

(955) 
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or 

ft x f 2 > 0 (95c) 

then f t and Pj t 2 is on the wrong side of the boundary, the previous 77 line. If 
all the points at the k=2 line pass this test, the second test should be run before 
checking points at the next k level. This first test compares each point in the 
grid with points having neighboring j and k values. Thus, it uses relatively little 
computational time per point, but each point must be checked. 

If the first test checks out the second test should be to check whether the k 
line intersects a previous k line. Figure 64 shows an example of this defect and the 
main reason to check for it. For strong convexities, when distances between the 
apex and its two adjoining points on the boundary are significantly different, the 
hyperbolic generator can cause the radial line emerging from the apex, the £ line 
(the radial line near the bottom in figure 64), to swing far enough over such that 
intersections occur. This can be checked by testing the 77 line segment at each 
point for intersection with the next lower 77 line. This check, like the previous 
one, though it checks each point in the grid, makes only one or two checks for 
intersection per point. 

A third breakdown, differing from the previous two in that there is no relation 
between the j values of the points compared, occurs when different parts of a 
boundary (the kmax line), or boundaries of different grids, intersect. The simplest 
procedure is to test each boundary line segment for intersection with all other 
boundary line segments. An improvement on this procedure is the following. Since 
intersection occurs only for lines whose x and y values overlap for some portion 
of the lines, computational time can be reduced, if the the check for intersection 
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is made only after the simpler checks for x and y overlap are made. This can be 
done by first sorting either x or y into a monotonic sequence and examining the 
new arrangement of the old j values. When successive j values differ by more than 
one, indicating overlap, the other overlap check can be made. In other words, if 
the original sequence is x(j), the new sequence can be denoted as x(j(i)). Then if 
j(i) and j(i+l) differ by more than one, there is overlap of x values of two separate 
line segments and the next test, to check for overlap of y values is made. If there 
is y overlap the check for intersection can be made. Because only points on the 
outer boundary are checked, the check for boundary intersection takes about the 

of time as the previous two procedures the check for 
order of time as the previous two. 


same order of magnitude 
breakdown is of the same 
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Grid Breakdown Due to Strong Convexity 
and Unequal Spacing Along Body 



Figure 64: tau=.5, cam=2.0, detail of hg5y20.as used as initial grid 
File : hg5y20.ptdps, bold line indicates body/boundary 
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8 Derivation of the Navier-Stokes Equations 


The governing equations in fluid dynamics are based on the three conservation 
laws of mass, momentum, and energy. These equations can be derived in the 
following way [10], [6] using a control- volume of fixed dimensions within which 
the three laws are analyzed. 


«*■+ 

A 



Figure 65: Fixed control volume showing area vector A and the three surface 
stresses at one of the two surfaces normal to the x axis. 

For the mass, or continuity, equation the change in mass per unit time, dm/dt 
within the control volume ( vol ), is balanced by the net sum of mass entering and 
leaving the volume in unit time. Since density, p, equals dm/dvol , the rate of 
change in mass within the volume is 

(dm\ _ d J dm 

\ dt / within dt 
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dfp dvol 


I 


dt 

dp 


dt 


dvol. 


(97) 

(98) 


The mass entering the volume per unit time can be represented as 


fdm\ 

\ dt / entering 


— <fc . pu • dA. 

/surface 


(99) 


The minus sign is required because A is directed away from the body. Com- 
bining equations (98) and (99) produces 


= ~ Lface^ 


( 100 ) 


— ♦ — f -» — * 

Using Gauss’ theorem, J G • dA = / V • Gdvol , equation (100) becomes 


J dvol = — J V • (p u)dvol , (101) 

which leads to the point form, since integration is the same for all terms, 

^ + V-(pu) = 0 (102) 

The first term on the left exists when mass density within the control volume 
is unsteady. The second term indicates the rate of mass per unit volume building 
up within the control volume due to differences in density and/or flow rates at 
the various surfaces. When the first term is zero, indicating no change in mass 
density at a point, the second term must also be zero. 

The momentum equation is a vector equation, containing three scalar equa- 
tions, one for each direction. It is derived from Newton’s second law: the rate of 
change in momentum for a body is equal to the forces acting on the body. For 
the fixed control volume, with mass flowing through it, the rate of change of mo- 
mentum (of the mass within it), is equal to the net rate momentum is convected 
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across its entire surface by the flow, plus the sum of the forces acting on the mass 
within it 

§ = § + * + p ‘ < 103 > 

where P{ is the momentum within the control volume, P e is the momentum 

— * — * 

convected across the surface, Fi represents body forces, such as gravity, and F, 

represents surface forces such as pressure and friction. 

Since dP = u dm, because u may vary within the control volume, and dm = 

pdvol, the rate of change of momentum within the volume is 

P i= [dPi= [ pudvol (104) 

J Jvol 

The momentum convected across the volume surface, per unit time, (covering 

a distance dr) is found by multiplying the differential of momentum at the sur- 

— # 

face with the dot product of the flow velocity u and the surface element A, and 
integrating over the whole surface 


dPc 

dt 


d(Jvoi pudvol) 

dt 

(105) 

d(S r Jarea pudr - dA) 

dt 

(106) 

r dr -* 

(107) 

_ / . dA 

Jarea H dt 

— / puu • dA 

Jarea 

(108) 

— I V • puudvol 

Jvol 

(109) 


The last change is by Gauss’ theorem. As before, the negative sign in front of 
the integral is due to the direction of the area element. 
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When these momentum terms axe on one side of the equation they then equal 
the sum of the forces acting on the mass in the control volume: body forces, such 
as gravity, and surface forces, or stresses, such as pressure and friction. 

For the control volume the body forces are integrated over the entire volume. 
If Fb is the total of body forces on the mass m in the volume then fo = Fb/m is 
the force per unit mass. And since density is defined as p = dm/dvol , pfb is the 
force per unit volume. Then the total force due to body forces is 

Fb = f d Pfbdvol. ( 110 ) 

Surfaces forces can be divided into two distint types: normal stresses and 
tangential stresses. Normal stresses, denoted as are pressure forces and 

act normal to the surface in question. The subscript indicates both the surface 
being acted upon and the direction of the force, which is normal to the surface. 
So B x acts on a surface area dA lying on a y-z plane (the line normal to a y-z 
plane lies in the x direction). 

Tangential forces are shearing, or friction, stresses, such as those associated 
with viscosity. At a surface there are two independent tangential directions, there- 
fore an indication of this must be made in denoting the symbol for this stress: 

A 

Tut = Tikk. The first subscript, like that for the normal stress, denotes the plane 
the surface lies on. The second subscript indicates the direction of the force. So 
the two tangential forces acting at the y-z surface dAx are r xy = r xy y acting in 
the y direction, and t xz = r xz z, acting in the z direction. 

The normal and shearing stresses form a stress tensor which is symmetric 
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) 


(in) 

( 112 ) 


The surfaces forces are integrated over the entire surface. Since the surface 

stresses form a tensor, its transformation by Gauss’ theorem produces the follow- 

— # 

ing equation for the sum of the surfaces forces F t acting on the control volume 


F a = f S-dA (113) 

Jarea 

= f V-Sdvol (114) 

where 

V-5 = ^ i,j = 1,2,3. (115) 

The momentum equation is then, after removing the integration signs, 

+ v ■ (pun) = V-S + pf b (116) 

and the x component is 

d(pUx) , ^7 ( d<j x dT xy 9 t xz j j ^ 

-^- + V-{pu,u) = — + — + — + PU (117) 

Since <t is the negative of pressure this equation becomes 


d(pu x ) 

dt 


+ V • ( pu x u ) - 


dPx drxy dr xz 

dx dy dz 


+ Pfb 


(118) 
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The equations for the y and z components of momentum are found in a similar 
way. / T xy T xz \ 

The energy equation derives from the first law of thermodynamics 

dE = dQ — dW (119) 

In the case of a control volume with mass passing through it, the rate Ei, the 
energy within the control volume, is changing, equals the rate energy, E c , is being 
convected in by the flow, plus the rate heat, Q conducted’ * s con( lucted in (ignoring 
heat radiated in for flows of sufficient speed), plus the rate work Wb is being done 
on the mass within the volume due to body forces, plus the rate work, W ap , is 
being done on the mass due to surface pressure forces, plus the rate work, W av , is 
being done on the mass due to surface viscous forces 


dEj dEc . ^Qconducted , dWb dWjp 

dt dt dt dt dt dt 


( 120 ) 


The energy per unit mass is p{e + l/2u 2 + e p ), with e, 1/2 u 2 and e p representing 
internal, kinetic, and potential energy respectively; u 2 = u • u. Similar to the 
previous equations for mass and momentum, the energy terms are replaced by 


= Lai P ( 6i + 1 / 2 ^ 2 + e Pi) dv ° l ( 121 ) 

^ = -I />(e + l/2n 2 + e p )w-di (122) 

dt Jarea 

= - f vol V • p(e c + l/2u 2 + epc)udvol. (123) 

By Fourier’s law, heat conduction per unit area is q— —kdT/dn, where k is a 
proportionality constant, T is temperature, n is the direction normal to the area, 
and the negative sign indicates that heat flux is from higher to lower temperatures. 


104 


In vector form, and integrating over the entire surface of the control volume the 
heat term in equation (120) becomes 


^conducted = _ I _ k ^ . ( 124 ) 

dt Jarea dn ' 

= <f k^n-dA ( 125 ) 

Jarea (jxi ' 

= L*< k f*>** ( 126 > 

= f vol V • (kV ■ T)dvol. ( 127 ) 

The work W b done by the body forces is found by integrating these forces over 
the entire volume. Since W = f dW — Fdx and F = / pfdvol, then pf ib is the 
body force per unit volume in the X{ direction and the energy due to this body 
force is 


dW b 

dt 


[dm 

dt 

I pfibdxidvol 


dt 

pfibdxidvol 

dt 

PfibdX{ 


I 
/ 

J pfxbUidvol. 


dt 


dvol 


This leads to 


dW b 

dt 


= J pfb • udvol 


( 128 ) 

( 129 ) 

( 130 ) 

( 131 ) 

( 132 ) 


( 133 ) 


Alternatively, this energy term (in fact all the energy terms associated with the 
momentum equation ) can be obtained from the momentum equation by taking 
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the dot product of the momentum change due to the body forces with the velocity 
vector-of’the flow ' '' ‘ ’-.- *•**—*> 


Fb'U = / ofb- udvol 

Jvol 


(134) 


Applied to the surface force term in the momentum equation, the technique 
produces 


dW, n dW, 


and that due to the x component of change is 

dpx . &T X y dT xz 

7\ P X H n H rj U X + 

ox oy oz 


(135) 


(136) 


The energy equation is then, after dropping the integration signs to get the 
differential form, and moving the convection energy term to the left, 


d[l( e 4 + ~2~ + e P»)l 

dt 


vf 


+ ^ * [p( e c + — + epc)]u = 

V • ( kVT ) + pf b • iT + (V • 5) • it. (137) 
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9 Results and Conclusions 


Steger’s original algorithm for solving the combined elliptic-hyperbolic equations, 
using fixed boundaries, was shown to work well on certain bodies, for certain 
initial grids. As concavity increased, the initial grids generated by the hyperbolic 
generator, hyg2d, broke down, necessitating the use of some other initial-grid 
generator. This latter generator, not being a function of the body, spaced points 
on the boundary irrespective of body shape, causing loss of orthogonality in some 
areas. 

The author then made modifications to the original algorithm to increase the 
range of bodies that the combined solver can handle. The first modification was 
to float the outer boundary. That is, at each iteration, after the interior of the grid 
was improved, the hyperbolic component, using the improved grid, generated a 
new and generally improved outer boundary. This allowed the use of initial grids 
with defects in the outer boundary. But, though an increased body complexity 
could now be dealt with, still, for strong enough concavities this algorithm breaks 
down. Additionally, as concavity increased, /z, the parameter used to control the 
amount of elliptic or hyperbolic, became more limited in the values it could take 
on. Only values within this range would allow for convergence to solution. Finding 
these values became more difficult as cavities grew in depth. 

A second modification allows the use of a different kind of initial grid, one 
with few or no defects. This is a grid of compressed grid cells where the length 
of cells in the rj direction is very small compared to the length in the £ direction. 
By starting with such a defect-free grid and allowing grid cells to grow with each 
iteration, combined with the floating outer boundary it is possible to generate a 
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grid about more complex (higher concavity) bodies. Additionally fi can be more 
arbitrarily picked, affording more variability between hyperbolic and elliptic. 

In the last modification to the algorithm, boundary movement within a cavity 
is modified. Instead of being generated by the hyperbolic alone it is modified to 
some degree by a forcing function that moves boundary points about a a chosen 
boundary point so that corresponding boundary points merge. This causes a cut 
to form within the cavity, made up of the outer boundary. This modification is 
not complete. There still remains the problem of forcing the point at which the 
outer boundary joins itself to move far enough from the body such that the entire 
relevant flow field is covered. 

When this is done a grid about a two-dimensional body with deep cavities 
may be achieved. Then, it should also be possible to generate a grid about mul- 
tiple bodies. By attaching them with imaginary lines, such as to make one inner 
boundary, the cavities so formed may be manageable. This idea is demonstrated 
in figures 66 through 69. Figure 66 shows two bodies connected by a line, which 
forms a continuous inner boundary that outlines the two bodies. Figure 67 is a 
grid generated about this new body. Figures 68 and 69 are detailed views of the 
area surrounding the connection line. The main point here is that there is little 
if any breakdown of the grid, no intersections are noticeable. 

Several methods of identifying grid breakdown have been described. Although 
the logic has been coded, and the first two techniques are in the grid generator 
program ehgrd, there remains yet the problem of using this information to control 
fi in order to eliminate the breakdown. One possibility is to simply turn the 
elliptic on (set fi high) wherever grid breakdown is detected. This, though, may 
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hamper using fi for resolution of flow variables. 
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Drawing One Grid About Two bodies by Connecting 
Them with a Line to Form One Body 



Figure 66: Airfoil and slat connected by a line 
file : f 1 4fp2.ps 



Figure 67: 100th pass of the combined solver, boundary is 
floated and grid cell volumes increase 5 percent per iteration 
File : fl 4f3.1 Ops 
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Drawing One Grid About Two bodies by Connecting 
Them with a Line to Form One Body 



Figure 68: 100th pass, detail of Final grid in Figure 66 
File : fMO.lOptaps 
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Drawing One Grid About Two bodies by Connecting 
Them with a Line to Form One Body 



Figure 69: 1 00th pass, detail of final grid in Figure 66 
File : fHO.lOptbps 
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10 Notation and Definitions 


u x = 

^XX = 


du 

dx 

d 2 u 

dx 2 


( 138 ) 

( 139 ) 



Linear Finite Difference Operators 


Operator 

Symbol 


Difference Representation 


Forward 

A* 

A x u 

= 

u{x + Ax) — tt(x) 

(140a) 

(Indexed) 


A x u 

= 

u(x i+ i) -u(xi) 

(140b) 




= 

U i+i - Ui 

(140c) 

Backward 

V* 

v*« 

— 

u(x) — u(x — A x ) 

(141a) 

(Indexed) 



= 

u(xi) - u(x,_i) 

(141b) 




= 

Ui - Ui- 1 

(141c) 

Central 

s x 

5 x u 

= 

u(x + Ax) — u(x — Ax) 

(142a) 

(Indexed) 


S x u 

= 

u(x i+1 ) - ti(x,_i) 

(142b) 




— 

«*+l - u i-l 

(142c) 

Central 

S x 

S x u 

— 

u(x + |Ax) — u(x — |Ax) 

(143) 

Central 

SI 

Slu 

— 

S x (S x u) 

(144a) 




= ■ 

u(x + Ax) — 2u(x) + u{x — Ax) 

(144b) 

(Indexed) 


S\u 

= 

tt(x i+ i) - 2u(xi) + u(xi_i) 

(144c) 




— 

Ui +1 - 2 Ui + Ui-l 

(144d) 


2nd order accurate, 1st partial derivative 


113 









= 2A X + ° ^ 

(145) 

. 6 x u 

u~ = 

* 2Ax 

(146) 

2nd order accurate, 2nd partial derivative 


f) 2 1! 

(147) 

6 2 u 

(148) 

• X 

- (Ax) 2 


iiu = u(x,y), and if only a finite number of x and y values are required, as is 
generally the case in numerical analysis, these values can be indexed as 

X = x(l) = X\ 

y = y(l) = yi 1 = 1,2, .. Amax. (149) 

u = u{l) = U\ 

Additionally, if the spacing in x and y values is constant a rectangular grid of 
points is formed, allowing double indexing of u 


X = 

x(m) 

— &m 

m = 1 , 2 , 

, . . . mmax 

y = 

y(n) 

= l In 

n — 1 , 2 , 

. . . nmax 

u = 

u(m,n) 

— 




( 150 ) 
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But if this is not the case a transformation of variables may be performed that 
allows such a rectangular grid to exist When x = x(£,t)) and y = £ and 77 


can be indexed as 

f = £(j) = £j 3 = 1,2,... jmax, (151) 

77 = r](k) = rjk k = 1, 2, . . . kmax, (152) 

This leads to 

X = x(j, k) = x jtk , (153) 

V = y(j,k) = y j}k (154) 

u = u(j,k) = Uj } k (155) 

If u is some derivative with respect to £ or tj the finite difference equivalent 
can be stated as 

imax 

v z(j->k) = £ cnv(j + i,k) (156) 

i—imin 

imax 

V v (. jyk)= £ biV { j,k + i) (157) 

i—imin 

Now if u happens to be a derivative with respect to x or y the chain rule can 
be used to transform into £,77 variables so that 

U x — -f- UijTlx (158) 
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Also by use of the metric transformation equations £ x and r) x can be changed to 

the u’s neighboring u(j,k), so that if v x = u, v x = v x (j, k), and, for a given value 
of j and k, j°, k° 

imax 

v x (j> k )- £ «<»(/ + *,fc) (159) 

where a,- is some rational constant and i may take on negative values. 


Convexity and Concavity 

Convexity refers to the convex portions of a body. Concavity refers to the 
relative depth and curvature of the concave portions of a body. The algorithm 
used to generate the two-dimensional bodies used in this paper controlled the 
body shape by a parameter named cam. This paramter when set to zero, as in 
figure 70, produces a body with no concave portion. As cam increases the bodies 
produced have what can be called stronger concavities as shown in figures 71, 72 
and 73. 
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The Relation of the Parameter Cam to Concavity 
' parameter tau = 0.5 



Figure 70: cam= 0.0 



Figure 72: cam= 0.80 



Figure 71: cam=0.40 



Figure 73 : cam= 1.0 
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